import os
import opensim as osim
import numpy as np
import pandas as pd

import matplotlib.pyplot as plt

import sys
sys.path.append(os.path.abspath('../../../analysis/'))

import importlib
import jam_analysis
from jam_analysis import JamAnalysis
from jam_analysis import get_h5_output

start_time = 1.16
start_pad = 0.0
stop_time = 2.36

GRF Data

def get_osim_table_results(table, return_pandas=True):
    data_array = table.getMatrix().to_numpy()
    colnames = table.getColumnLabels()
    time = np.asarray(table.getIndependentColumn())
    
    if return_pandas is True:
        df = pd.DataFrame()
        df['time'] = time
        for idx, col_name in enumerate(colnames):
            df[col_name] = data_array[:, idx]
        return df
    
    else:
        return {
            'time': time,
            'colname': colnames,
            'data': data_array
        }
    

grf_mot_file = '../../../../models/knee_healthy/experimental_data/motion_analysis/overground_17_grf.mot'

grf_data = osim.TimeSeriesTable(grf_mot_file)
df_grf = get_osim_table_results(grf_data)
# grf_array = mgrf_data.getMatrix().to_numpy()
# grf_colnames = grf_data.getColumnLabels()
# time = np.asarray(grf_data.getIndependentColumn())



start_idx = np.where(df_grf['time'] == start_time)[0][0]
end_idx = np.where(df_grf['time'] == stop_time)[0][0]

grf_freq = 2000
grf_time = np.arange(start_time, stop_time, 1/grf_freq)


forceplates = (1, 2, 3)
parameters = (('v', 'force'), ('p', 'force'), ('', 'torque'))
axes = ('x', 'y', 'z')
data_string = 'ground_{parameter_name}_{plate_number}_{parameter}{axis}'

fig, ax = plt.subplots(3, 1, figsize=(16,18))

for idx, (param_num, param_name) in enumerate(parameters):
    for fp in forceplates:
        for axis in axes:
            colname = data_string.format(parameter_name=param_name, plate_number=fp, parameter=param_num, axis=axis)
#             col_idx, colname_ = [(idx, x) for idx, x in enumerate(grf_data['colname']) if x == colname][0]
            
            ax[idx].plot(grf_time, df_grf[colname][start_idx:end_idx], label=colname)
            ax[idx].set_title(colname)
            
ax[0].set_title('Ground Reaction Force (N)')
ax[0].legend()
ax[1].set_title('Ground Reaction CofP (m?)')
ax[1].legend()
ax[2].set_title('Ground Reaction Torque (Nm)')
ax[2].legend()
                            

<matplotlib.legend.Legend at 0x7fb4d52b5190>

COMAK Inverse Kinematics Results

secondary_constraint_function_file = './results/comak-inverse-kinematics/secondary_coordinate_constraint_functions.xml'

ik_result_dir = './results/comak-inverse-kinematics'
results_basename = "walking" # for file sharing

func_set = osim.FunctionSet(secondary_constraint_function_file)

df_sweep_results = get_osim_table_results(
    osim.TimeSeriesTable(os.path.join(
        ik_result_dir, 
        results_basename + '_secondary_constraint_sweep_states.sto')
                        )
)

# sweep_cols = [
#     '/jointset/knee_r/knee_add_r/value',
#     '/jointset/knee_r/knee_rot_r/value',
#     '/jointset/knee_r/knee_tx_r/value',
#     '/jointset/knee_r/knee_ty_r/value',
#     '/jointset/knee_r/knee_tz_r/value',
#     '/jointset/pf_r/pf_flex_r/value',
#     '/jointset/pf_r/pf_rot_r/value',
#     '/jointset/pf_r/pf_tilt_r/value',
#     '/jointset/pf_r/pf_tx_r/value',
#     '/jointset/pf_r/pf_ty_r/value',
#     '/jointset/pf_r/pf_tz_r/value',
# ]

sweep_flex = np.rad2deg(df_sweep_results['/jointset/knee_r/knee_flex_r/value'])

/miniconda/lib/python3.7/site-packages/ipykernel_launcher.py:10: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead.  To get a de-fragmented frame, use `newframe = frame.copy()`
  # Remove the CWD from sys.path while we load stuff.

fig, ax = plt.subplots(4, 3, figsize=(16,18))

for i in range(func_set.getSize()):
    func = osim.GCVSpline.safeDownCast(func_set.get(i))
    name = func.getName()
    
    X = np.zeros(func.getNumberOfPoints())
    Y = np.zeros(func.getNumberOfPoints())
    
    for pt_idx in range(func.getNumberOfPoints()):
        X[pt_idx] = func.getX(pt_idx)
        Y[pt_idx] = func.getY(pt_idx)
    
    if ('tx' in name) or ('ty' in name) or ('tz' in name):
        ax[i%4, i//4].plot(sweep_flex, df_sweep_results[name + '/value'], label='Sweep Simulation')
        ax[i%4, i//4].plot(np.rad2deg(X), Y, label='Constraint Function')
        ax[i%4, i//4].set_ylabel('Translation [mm]')        
    
    else:
        ax[i%4, i//4].plot(sweep_flex, np.rad2deg(df_sweep_results[name + '/value']), label='Sweep Simulation')
        ax[i%4, i//4].plot(np.rad2deg(X), np.rad2deg(Y), label='Constraint Function')
        ax[i%4, i//4].set_ylabel('Angle [^o]')
    ax[i%4, i//4].set_title(name)
    ax[i%4, i//4].set_xlabel('Knee Flexion [^o]')
    ax[i%4, i//4].legend()

plt.tight_layout()

COMAK Results

path_comak_h5_files = [f'./results/joint-mechanics/{results_basename}.h5',]

fontsize=20
linewidth=5

jam = JamAnalysis()
jam.jam_analysis(path_comak_h5_files)

Plot Secondary Kinematics

sec_coords = [
    'knee_flex_r','knee_add_r','knee_rot_r',
    'knee_tx_r','knee_ty_r','knee_tz_r',
]

rows = 2
cols = 3
fig, ax = plt.subplots(rows, cols, figsize=(16*3, 4*12))
for idx, outcome in enumerate(sec_coords):
    row_ = idx // cols
    col_ = idx % cols
    for file_idx in range(len(path_comak_h5_files)):
        ax[row_, col_].plot(jam.coordinateset[outcome]['value'][:,file_idx], linewidth=linewidth, label=file_idx)
    ax[row_, col_].set_title(outcome, fontsize=fontsize*3)
    ax[row_, col_].xaxis.set_tick_params(labelsize=fontsize*2)
    ax[row_, col_].yaxis.set_tick_params(labelsize=fontsize*2)
    if idx == 0:
        ax[row_, col_].legend(fontsize=fontsize*2)

fig.suptitle('Tibiofemoral Secondary Kinematics', fontsize=fontsize*4)
fig.tight_layout(rect=[0, 0.03, 1, 0.95])

sec_coords = [
    'pf_flex_r','pf_rot_r','pf_tilt_r',
    'pf_tx_r','pf_ty_r','pf_tz_r'
]

rows = 2
cols = 3
fig, ax = plt.subplots(rows, cols, figsize=(16*3, 4*12))
for idx, outcome in enumerate(sec_coords):
    row_ = idx // cols
    col_ = idx % cols
    for file_idx in range(len(path_comak_h5_files)):
        ax[row_, col_].plot(jam.coordinateset[outcome]['value'][:,file_idx], linewidth=linewidth, label=file_idx)
    ax[row_, col_].set_title(outcome, fontsize=fontsize*3)
    ax[row_, col_].xaxis.set_tick_params(labelsize=fontsize*2)
    ax[row_, col_].yaxis.set_tick_params(labelsize=fontsize*2)
    if idx == 0:
        ax[row_, col_].legend(fontsize=fontsize*2)

fig.suptitle('Patellofemoral Secondary Kinematics', fontsize=fontsize*4)
fig.tight_layout(rect=[0, 0.03, 1, 0.95])

Muslce Activations

list_all_muscles = list(jam.forceset['Muscle'].keys())

muscle_outcomes = [
    'activation',
#     'actuation',
#     'active_fiber_force',
#     'passive_fiber_force',
#     'tendon_strain',
#     'normalized_fiber_length',
#     'length'
]


cols = 3
rows = len(list_all_muscles) // cols + 1

fig, ax = plt.subplots(rows, cols, figsize=(16*cols, (rows+1)*12))

fig.suptitle('Muscle Activation Outputs - ALL MUSCLES', fontsize=fontsize*3.5)

muscle_outcome = muscle_outcomes[0]
# ACTUATION
for muscle_idx, muscle in enumerate(list_all_muscles):
    row_ = muscle_idx // cols
    col_ = muscle_idx % cols
#     for idx, muscle in enumerate(muscles):
    jam.plot_muscle_output(
        muscle, 
        muscle_outcome, 
        ax=ax[row_, col_], 
        linewidth=5, 
        fontsize=fontsize*3,
        label=muscle
    )
    ax[row_, col_].set_ylabel('Activation', fontsize=fontsize*3)
    ax[row_, col_].xaxis.set_tick_params(labelsize=fontsize*2)
    ax[row_, col_].yaxis.set_tick_params(labelsize=fontsize*2)
    ax[row_, col_].set_title(muscle, fontsize=fontsize*3)
    if muscle_idx == 2:
        ax[row_, col_].legend(fontsize=fontsize*2)
        
fig.tight_layout(rect=[0, 0.03, 1, 0.95])

knee_muscle_names = [
    'vasint_r','vaslat_r','vasmed_r','recfem_r',
    'bflh_r','bfsh_r','semimem_r','semiten_r',
    'gasmed_r','gaslat_r'
]


muscle_outcomes = [
    'activation',
#     'actuation',
#     'active_fiber_force',
#     'passive_fiber_force',
#     'tendon_strain',
#     'normalized_fiber_length',
#     'length'
]

cols = 3
rows = len(knee_muscle_names) // cols + 1

fig, ax = plt.subplots(rows, cols, figsize=(16*cols, (rows+1)*12))

fig.suptitle('Muscle Activation Outputs - KNEE MUSCLES', fontsize=fontsize*3.5)

muscle_outcome = muscle_outcomes[0]
# ACTUATION
for muscle_idx, muscle in enumerate(knee_muscle_names):
    row_ = muscle_idx // cols
    col_ = muscle_idx % cols
#     for idx, muscle in enumerate(muscles):
    jam.plot_muscle_output(
        muscle, 
        muscle_outcome, 
        ax=ax[row_, col_], 
        linewidth=5, 
        fontsize=fontsize*3,
        label=muscle
    )
    ax[row_, col_].set_ylabel(muscle_outcome, fontsize=fontsize*3)
    ax[row_, col_].xaxis.set_tick_params(labelsize=fontsize*2)
    ax[row_, col_].yaxis.set_tick_params(labelsize=fontsize*2)
    ax[row_, col_].set_title(muscle, fontsize=fontsize*3)
    if muscle_idx == 2:
        ax[row_, col_].legend(fontsize=fontsize*2)
        
fig.tight_layout(rect=[0, 0.03, 1, 0.95])

MUSCLE FORCES

knee_muscle_names = [
    'vasint_r','vaslat_r','vasmed_r','recfem_r',
    'bflh_r','bfsh_r','semimem_r','semiten_r',
    'gasmed_r','gaslat_r'
]


muscle_outcomes = [
#     'activation',
    'actuation',
#     'active_fiber_force',
#     'passive_fiber_force',
#     'tendon_strain',
#     'normalized_fiber_length',
#     'length'
]

cols = 3
rows = len(knee_muscle_names) // cols + 1

fig, ax = plt.subplots(rows, cols, figsize=(16*cols, (rows+1)*12))

fig.suptitle('Muscle Forces - KNEE MUSCLES', fontsize=fontsize*3.5)

muscle_outcome = muscle_outcomes[0]
# ACTUATION
for muscle_idx, muscle in enumerate(knee_muscle_names):
    row_ = muscle_idx // cols
    col_ = muscle_idx % cols
#     for idx, muscle in enumerate(muscles):
    jam.plot_muscle_output(
        muscle, 
        muscle_outcome, 
        ax=ax[row_, col_], 
        linewidth=5, 
        fontsize=fontsize*3,
        label=muscle
    )
    ax[row_, col_].set_ylabel('Force (N)', fontsize=fontsize*3)
    ax[row_, col_].xaxis.set_tick_params(labelsize=fontsize*2)
    ax[row_, col_].yaxis.set_tick_params(labelsize=fontsize*2)
    ax[row_, col_].set_title(muscle, fontsize=fontsize*3)
    if muscle_idx == 2:
        ax[row_, col_].legend(fontsize=fontsize*2)
        
fig.tight_layout(rect=[0, 0.03, 1, 0.95])

Ligament Forces

# lig_base_names = ['ACL','PCL','MCL','LCL','PT','ITB','pCAP','mPFL','lPFL']
# list_ = []
# for lig in jam.forceset['Blankevoort1991Ligament'].keys():
#     for lig_base_name in lig_base_names:
#         if lig_base_name in lig:
#             list_.append(lig)

# [x for x in jam.forceset['Blankevoort1991Ligament'].keys() if ]

# Ligament Forces
# ligament_names = [
#     'MCLd','MCLs',
#     'LCL',
#     'ACLam', 'ACLpl',
#     'PCLal', 'PCLpm',
#     'PT',
#     'ITB',
#     'pCAP'
# ]

# lig_base_names = {'ACL','PCL','MCL','LCL','PT','ITB','pCAP','mPFL','lPFL'};

lig_base_names = ['ACL','PCL','MCL','LCL','PT','ITB','pCAP','mPFL','lPFL']
ligament_names = [x for x in jam.forceset['Blankevoort1991Ligament'].keys() for y in lig_base_names if y in x]

cols = 3
rows = len(ligament_names) // cols + 1


fig, ax = plt.subplots(rows, cols, figsize=(16*cols, (rows+1)*12))
fig.suptitle('Ligament Forces', fontsize=fontsize*4)

for idx, outcome in enumerate(ligament_names):
    row_ = idx // cols
    col_ = idx % cols
    
    fibers = [x for x in jam.forceset['Blankevoort1991Ligament'].keys() if outcome in x]
    
    
    for file_idx in range(len(path_comak_h5_files)):
        data = np.zeros(jam.forceset['Blankevoort1991Ligament'][fibers[0]]['total_force'].shape[0])
        for fiber in fibers:
            data += np.squeeze(jam.forceset['Blankevoort1991Ligament'][fiber]['total_force'][:, file_idx])
        ax[row_, col_].plot(data, linewidth=5, label=file_idx)
    ax[row_, col_].set_title(outcome, fontsize=fontsize*3)
    ax[row_, col_].xaxis.set_tick_params(labelsize=fontsize*2)
    ax[row_, col_].yaxis.set_tick_params(labelsize=fontsize*2)
    if idx == 0:
        ax[row_, col_].legend(fontsize=fontsize*2)
    
fig.tight_layout(rect=[0, 0.03, 1, 0.95])

Ligament Strain

# Ligament Forces
# ligament_names = [
#     'MCLd','MCLs',
#     'LCL',
#     'ACLam', 'ACLpl',
#     'PCLal', 'PCLpm',
#     'PT',
#     'ITB',
#     'pCAP'
# ]

# lig_base_names = {'ACL','PCL','MCL','LCL','PT','ITB','pCAP','mPFL','lPFL'};

lig_base_names = ['ACL','PCL','MCL','LCL','PT','ITB','pCAP','mPFL','lPFL']
ligament_names = [x for x in jam.forceset['Blankevoort1991Ligament'].keys() for y in lig_base_names if y in x]

cols = 3
rows = len(ligament_names) // cols + 1


fig, ax = plt.subplots(rows, cols, figsize=(16*cols, (rows+1)*12))
fig.suptitle('Ligament Strains', fontsize=fontsize*4)

for idx, outcome in enumerate(ligament_names):
    row_ = idx // cols
    col_ = idx % cols
    
    fibers = [x for x in jam.forceset['Blankevoort1991Ligament'].keys() if outcome in x]
    
    
    for file_idx in range(len(path_comak_h5_files)):
        data = np.zeros(jam.forceset['Blankevoort1991Ligament'][fibers[0]]['strain'].shape[0])
        for fiber in fibers:
            data += np.squeeze(jam.forceset['Blankevoort1991Ligament'][fiber]['strain'][:, file_idx])
        ax[row_, col_].plot(data, linewidth=5, label=file_idx)
    ax[row_, col_].set_title(outcome, fontsize=fontsize*3)
    ax[row_, col_].xaxis.set_tick_params(labelsize=fontsize*2)
    ax[row_, col_].yaxis.set_tick_params(labelsize=fontsize*2)
    if idx == 0:
        ax[row_, col_].legend(fontsize=fontsize*2)
    
fig.tight_layout(rect=[0, 0.03, 1, 0.95])

Knee Contact Forces

cols = 3
rows = 5


fig, ax = plt.subplots(rows, cols, figsize=(16*cols, (rows+1)*12))
fig.suptitle('Contact Forces (N)', fontsize=fontsize*4)

for axis_idx, axis in enumerate(['X', 'Y', 'Z']):
    # TF Total
    ax[0, axis_idx].plot(
        jam.time,
        jam.forceset['Smith2018ArticularContactForce']['tf_contact']['tibia_cartilage']['total_contact_force'][:,axis_idx],
        linewidth=5
    )
    ax[0, axis_idx].set_title('TF Total Force ' + axis, fontsize=fontsize*3)
    # TF Medial
    ax[1, axis_idx].plot(
        jam.time,
        jam.forceset['Smith2018ArticularContactForce']['tf_contact']['tibia_cartilage'][4]['regional_contact_force'][:,axis_idx],
        linewidth=5
    )
    ax[1, axis_idx].set_title('Medial Tibia Force ' + axis, fontsize=fontsize*3)
    
    # TF Lateral
    ax[2, axis_idx].plot(
        jam.time,
        jam.forceset['Smith2018ArticularContactForce']['tf_contact']['tibia_cartilage'][5]['regional_contact_force'][:,axis_idx],
        linewidth=5
    )
    ax[2, axis_idx].set_title('Lateral Tibia Force ' + axis, fontsize=fontsize*3)
    
#     plot(jam.time,jam.forceset.Smith2018ArticularContactForce.pf_contact.patella_cartilage.total_contact_force(:,i))
    
    # PF Total
    ax[3, axis_idx].plot(
        jam.time,
        jam.forceset['Smith2018ArticularContactForce']['pf_contact']['patella_cartilage']['total_contact_force'][:,axis_idx],
        linewidth=5
    )
    ax[3, axis_idx].set_title('Patellofemoral Force ' + axis, fontsize=fontsize*3)
    
    # TF - Total, Medial, Lateral
    ax[4, axis_idx].set_title('All Tibiofemoral Forces ' + axis, fontsize=fontsize*3)
    ax[4, axis_idx].plot(
        jam.time,
        jam.forceset['Smith2018ArticularContactForce']['tf_contact']['tibia_cartilage']['total_contact_force'][:,axis_idx],
        linewidth=5,
        label='Total'
    )
    
    ax[4, axis_idx].plot(
        jam.time,
        jam.forceset['Smith2018ArticularContactForce']['tf_contact']['tibia_cartilage'][4]['regional_contact_force'][:,axis_idx],
        linewidth=5,
        label='Medial'
    )
    ax[4, axis_idx].plot(
        jam.time,
        jam.forceset['Smith2018ArticularContactForce']['tf_contact']['tibia_cartilage'][5]['regional_contact_force'][:,axis_idx],
        linewidth=5,
        label='Lateral'
    )
    
    total_ = (
        jam.forceset['Smith2018ArticularContactForce']['tf_contact']['tibia_cartilage'][4]['regional_contact_force'][:,axis_idx] 
        + jam.forceset['Smith2018ArticularContactForce']['tf_contact']['tibia_cartilage'][5]['regional_contact_force'][:,axis_idx] 
    )
    
    ax[4, axis_idx].plot(
        jam.time,
        total_,
        linewidth=5,
        linestyle='--',
        label='Total (Lat + Med)'
    )
    ax[4, axis_idx].legend(fontsize=40)

    
    
    for row_ in range(rows):
#         ax[row_, axis_idx].set_title(axis, fontsize=fontsize*3)
        ax[row_, axis_idx].xaxis.set_tick_params(labelsize=fontsize*2)
        ax[row_, axis_idx].yaxis.set_tick_params(labelsize=fontsize*2)
        

# for row_ in range(rows):
#     for col_ in enumerage(('X', 'Y', 'Z'))
#     ax[row_, col_].set_title(outcome, fontsize=fontsize*3)
#     ax[row_, col_].xaxis.set_tick_params(labelsize=fontsize*2)
#     ax[row_, col_].yaxis.set_tick_params(labelsize=fontsize*2)
    
    
    

Knee Contact Pressures

cols = 4
rows = 1


fig, ax = plt.subplots(rows, cols, figsize=(16*cols, (rows+1)*9))
fig.suptitle('Contact Presures (KPa??)', fontsize=fontsize*4)

    # TF Total
#     plot(jam.time,jam.forceset.Smith2018ArticularContactForce.tf_contact.tibia_cartilage.total_max_pressure)
#     plot(jam.time,jam.forceset.Smith2018ArticularContactForce.tf_contact.tibia_cartilage.total_mean_pressure)
#     title(['TF Total' ])
    
ax[0].plot(
    jam.time,
    jam.forceset['Smith2018ArticularContactForce']['tf_contact']['tibia_cartilage']['total_max_pressure'],
    linewidth=5,
    label='Max Pressure'
)
ax[0].plot(
    jam.time,
    jam.forceset['Smith2018ArticularContactForce']['tf_contact']['tibia_cartilage']['total_mean_pressure'],
    linewidth=5,
    label='Mean Pressure'
)
ax[0].legend(fontsize=40)
ax[0].set_title('Tibial Cartilage Total', fontsize=50)

ax[1].plot(
    jam.time,
    jam.forceset['Smith2018ArticularContactForce']['tf_contact']['tibia_cartilage'][4]['regional_max_pressure'],
    linewidth=5,
    label='Max Pressure'
)
ax[1].plot(
    jam.time,
    jam.forceset['Smith2018ArticularContactForce']['tf_contact']['tibia_cartilage'][4]['regional_mean_pressure'],
    linewidth=5,
    label='Mean Pressure'
)
ax[1].legend(fontsize=40)
ax[1].set_title('Medial Tibial Cartilage', fontsize=50)

ax[2].plot(
    jam.time,
    jam.forceset['Smith2018ArticularContactForce']['tf_contact']['tibia_cartilage'][5]['regional_max_pressure'],
    linewidth=5,
    label='Max Pressure'
)
ax[2].plot(
    jam.time,
    jam.forceset['Smith2018ArticularContactForce']['tf_contact']['tibia_cartilage'][5]['regional_mean_pressure'],
    linewidth=5,
    label='Mean Pressure'
)
ax[2].legend(fontsize=40)
ax[2].set_title('Lateral Tibial Cartilage', fontsize=50)

ax[3].plot(
    jam.time,
    jam.forceset['Smith2018ArticularContactForce']['pf_contact']['patella_cartilage']['total_max_pressure'],
    linewidth=5,
    label='Max Pressure'
)
ax[3].plot(
    jam.time,
    jam.forceset['Smith2018ArticularContactForce']['pf_contact']['patella_cartilage']['total_mean_pressure'],
    linewidth=5,
    label='Mean Pressure'
)
ax[3].legend(fontsize=40)
ax[3].set_title('Tibial Cartilage Total', fontsize=50)

for idx in range(4):
    ax[idx].xaxis.set_tick_params(labelsize=fontsize*2)
    ax[idx].yaxis.set_tick_params(labelsize=fontsize*2)
#     ax[0, axis_idx].set_title('TF Total Force ' + axis, fontsize=fontsize*3)

Knee COP

cols = 3
rows = 1

fig, ax = plt.subplots(rows, cols, figsize=(16*cols, (rows+1)*9))
fig.suptitle('Center of Pressure (mm)', fontsize=fontsize*4)

tf_med_cop = jam.forceset['Smith2018ArticularContactForce']['tf_contact']['tibia_cartilage'][4]['regional_center_of_pressure'] * 1000
ax[0].plot(tf_med_cop[:, 2], tf_med_cop[:, 0], linewidth=5)
ax[0].set_title('Medial Tibia', fontsize=50)
ax[0].set_xlabel('Lateral [mm]', fontsize=40)
ax[0].set_ylabel('Anterior [mm]', fontsize=40)
ax[0].set_xlim(-30, 0)
ax[0].set_ylim(-15, 15)

tf_lat_cop = jam.forceset['Smith2018ArticularContactForce']['tf_contact']['tibia_cartilage'][5]['regional_center_of_pressure'] * 1000
ax[1].plot(tf_lat_cop[:, 2], tf_lat_cop[:, 0], linewidth=5)
ax[1].set_title('Lateral Tibia', fontsize=50)
ax[1].set_xlabel('Lateral [mm]', fontsize=40)
ax[1].set_ylabel('Anterior [mm]', fontsize=40)
ax[1].set_xlim(0, 30)
ax[1].set_ylim(-15, 15)

pf_cop = jam.forceset['Smith2018ArticularContactForce']['pf_contact']['patella_cartilage']['total_center_of_pressure'] * 1000
ax[2].plot(pf_cop[:, 2], pf_cop[:, 1], linewidth=5)
ax[2].set_title('Patella', fontsize=50)
ax[2].set_xlabel('Lateral [mm]', fontsize=40)
ax[2].set_ylabel('Superior [mm]', fontsize=40)
ax[2].set_xlim(-15, 15)
ax[2].set_ylim(-15, 15)

for idx in range(3):
    ax[idx].xaxis.set_tick_params(labelsize=fontsize*2)
    ax[idx].yaxis.set_tick_params(labelsize=fontsize*2)

Inverse Dynamics

id_coords = [
    'hip_flex_r','hip_add_r','hip_rot_r',
    'knee_flex_r','knee_add_r','knee_rot_r',
    'knee_tx_r','knee_ty_r','knee_tz_r',
    'pf_flex_r','pf_rot_r','pf_tilt_r',
    'pf_tx_r','pf_ty_r','pf_tz_r',
    'ankle_flex_r','subt_angle_r','mtp_angle_r'
]

id_sto_file = './results/comak-inverse-dynamics/walking_inverse-dynamics.sto'

id_data = osim.TimeSeriesTable(id_sto_file)
df_id = get_osim_table_results(id_data)

cols = 3
rows = 6

fig, ax = plt.subplots(rows, cols, figsize=(16*cols, (rows+1)*9))
# fig.suptitle('Inverse Dynamics Results', fontsize=fontsize*4)

for idx, id_coord in enumerate(id_coords):
    row_ = idx // cols
    col_ = idx % cols
    
    if ('tx' in id_coord) or ('ty' in id_coord) or ('tz' in id_coord):
        type_ = 'force'
        unit = '[N]'
    else:
        type_ = 'moment'
        unit = '[Nm]'
    data_name = f'{id_coord}_{type_}'
    ax[row_, col_].plot(df_id['time'], df_id[data_name], linewidth=5)
    ax[row_, col_].set_title(data_name, fontsize=fontsize*3)
    ax[row_, col_].xaxis.set_tick_params(labelsize=fontsize*2)
    ax[row_, col_].yaxis.set_tick_params(labelsize=fontsize*2)
    ax[row_, col_].set_xlabel('Time [s]', fontsize=40)
    ax[row_, col_].set_ylabel(f'{type} {unit}', fontsize=40)

suptitle = plt.suptitle('Inverse Dynamics Results', fontsize=fontsize*4, y=1.02)

plt.tight_layout()